#Kris-Stella Trump kstrump@gmail.com 
#This is a data analysis file for experiment carried out July 2-9 2013 - using the inescapability paragraph to manipulate system justification and measure whether this affects adjustment to information regarding high income inequality.

library(foreign)
alldata <- read.csv("Exp4 - cleandata.csv", na.strings="")
names(alldata)


#Remove the 8 respondents marked as suspicious for giving nonsensical answers 
data <- alldata[is.na(alldata$suspicious),]

#Create the variables needed for analysis

#Perceptions

data$percmin <- apply(data[,c("est_1_1_TEXT", "est_2_1_TEXT", "est_3_1_TEXT", "est_4_1_TEXT", "est_5_1_TEXT", "est_6_1_TEXT")], 1, min)
data$percmax <- apply(data[,c("est_1_1_TEXT", "est_2_1_TEXT", "est_3_1_TEXT", "est_4_1_TEXT", "est_5_1_TEXT", "est_6_1_TEXT")], 1, max)
data$percratio <- (data$percmax/data$percmin)
data$percindex <- log(data$percratio)
summary(data$percratio)


#Preferences total
data$prefmin <- apply(data[,c("tot.pref_x1", "tot.pref_x2", "tot.pref_x3", "tot.pref_x4", "tot.pref_x5", "tot.pref_x6")], 1, min)
data$prefmax <- apply(data[,c("tot.pref_x1", "tot.pref_x2", "tot.pref_x3", "tot.pref_x4", "tot.pref_x5", "tot.pref_x6")], 1, max)
data$prefratio <- (data$prefmax/data$prefmin)
data$prefindex <- log(data$prefratio)
summary(data$prefratio)

#Check for outliers
plot(data$prefratio) #non-logged
plot(data$prefindex) #logged

#Partisanship with leaners

#to avoid problems with na, make them zero in party.ind only
data$party.ind[is.na(data$party.ind)] <- 0
data$republican <- c()
data$republican <- ifelse((data$party>3 | data$party.ind==3), 1, 0)
data$democrat <- c()
data$democrat <- ifelse((data$party<3 | data$party.ind==1), 1, 0)
#both indicators include complete independents in the omitted group


#BJW indicator
data$bjw <- (data$bjw_1 + data$bjw_2 + data$bjw_3 + data$bjw_4 + data$bjw_5 + data$bjw_6 + data$bjw_7 + data$bjw_8) / 8
summary(data$bjw)
sum(data$bjw==3.75, na.rm=T)
plot(density(data$bjw, na.rm=T))

#Treatment indicators

#Information
data$info <- (is.na(data$tr.table))==FALSE
summary(data$info)

#High/low escapability: treatment variable
data$threat <- (is.na(data$high.ine.t_3))==FALSE 
summary(data$threat)



#Descriptive statistics

summary(data$born)
#mean 30, median 27

summary(data$percratio)
#median 30, mean 1357
#three-four outliers. look at log transformed variable. 
summary(data$percindex)
exp(3.4)
exp(3.6)
# median 30, mean 37

sum(data$gender==1, na.rm=T)
219 / 589 #=0.37

data$degree <- data$educ>=4
sum(data$degree)
316/589

sum(data$democrat==1, na.rm=T)
361/589 #=0.61
sum(data$republican==1, na.rm=T)
117/589 #=0.20

sum(data$race==1, na.rm=T)
417/589 #=71% white
sum(data$race==2, na.rm=T)
39/589 #=0.07 black
sum(data$race==4, na.rm=T)
76/589 #=13% asian
data$white <- data$race==1


##### Balance table

balance <- zelig(threat ~ born + gender + white + republican +degree + bjw + percindex, model="ls", data=data)
summary(balance)

#Age (born+18 for age)
mean(data[data$threat==0,"born"], na.rm=T)
mean(data[data$threat==1,"born"], na.rm=T) 
#Gender
sum(data[data$threat==0 & data$gender==1,"gender"], na.rm=T)/sum(data$threat==0)
sum(data[data$threat==1 & data$gender==1,"gender"], na.rm=T)/sum(data$threat==1)
#Race
sum(data[data$threat==0 & data$white==1,"white"], na.rm=T)/sum(data$threat==0)
sum(data[data$threat==1 & data$white==1,"white"], na.rm=T)/sum(data$threat==1)
#BJW
mean(data[data$threat==0,"bjw"], na.rm=T)
mean(data[data$threat==1,"bjw"], na.rm=T) 
#Republican
sum(data[data$threat==0 & data$republican==1,"republican"], na.rm=T)/sum(data$threat==0)
sum(data[data$threat==1 & data$republican==1,"republican"], na.rm=T)/sum(data$threat==1)
#College degree
sum(data[data$threat==0 & data$degree==1,"degree"], na.rm=T)/sum(data$threat==0)
sum(data[data$threat==1 & data$degree==1,"degree"], na.rm=T)/sum(data$threat==1)
#Perception of inequality
mean(data[data$threat==0,"percindex"], na.rm=T)
mean(data[data$threat==1,"percindex"], na.rm=T)



#Treatment effects, Predicted values

#install.packages("Zelig")
library(Zelig)

#Show that threat impacts preferences
z.out1 <- zelig(prefindex ~ info + threat, model="ls", data=data)
summary(z.out1)

#Threat interacts with info: most of the raw threat effect occurs among those who received info.
z.out2 <- zelig(prefindex ~ info*threat, model="ls", data=data)
summary(z.out2)


#Threat, bjw and political orientation all influence desired inequality
z.out3 <- zelig(prefindex ~ info + threat + bjw + republican + percindex, model="ls", data=data)
summary(z.out3)

#Adding interaction
z.out4 <- zelig(prefindex ~ info*threat + bjw + republican + percindex, model="ls", data=data)
summary(z.out4)



#Check impact of threat on prior beliefs about inequality
z.out.est <- zelig(percindex ~ info + threat + bjw + republican, model="ls", data=data)
summary(z.out.est)
#none -> successful randomization



################Graphing the results

library(Hmisc)

###Graph: threat impact concentrated among those who received info

#Create expected values. 
x.out <- setx(z.out2, threat=c(0,1), info=c(0,1))
s.out <- sim(z.out2, x=x.out)
summary(s.out)
#Graph it

#order: 
#no info no threat, no info and threat, no threat and info, info and threat

mean <- c(2.14, 2.13, 2.57,3.02) 
lci <- c(1.94, 1.92, 2.37, 2.82)
uci <- c(2.33, 2.34,2.78, 3.21)
at <- c(3,7)

bjwresultslog <- as.data.frame(cbind(mean, uci, lci))
bjwresults <- exp(bjwresultslog)
bjwresults <- cbind(bjwresults,at)
mp <- barplot(bjwresults$mean,ann=F, ylim=c(0,26), xaxt="n", space=c(0.5,0,0.5,0),col=c("grey40","lightgrey"))


pdf("threatresults.pdf")
barplot(bjwresults$mean,ann=F, ylim=c(0,30), xaxt="n", space=c(0.5,0,0.5,0),col=c("grey40","lightgrey"))
axis(1,at=c(0,10),lwd.ticks=0, labels=c("",""))
errbar(mp,bjwresults$mean,bjwresults$uci,bjwresults$lci,add=T, pch=NA)
legend("topleft", legend=c("Control", "System justification"), fill=c("grey40","lightgrey"))
title(main="Recommended income ratios", ylab="Income ratio")
mtext(c("No information", "Information condition"), side=1, at=c(1.5,4), line=0)
dev.off()



#Should government reduce income differences? 
z.out.red <- zelig(opinion_3 ~ info + threat + republican + bjw, model="ls", data=data)
summary(z.out.red)
z.out.redx <- zelig(opinion_3 ~ info*threat + republican + bjw, model="ls", data=data)
summary(z.out.redx)
#no impact of either treatment

#Income differences are too large
z.out.diff <- zelig(opinion_1 ~ info + threat + republican + bjw, model="ls", data=data)
summary(z.out.diff)
z.out.diffx <- zelig(opinion_1 ~ info*threat + republican + bjw, model="ls", data=data)
summary(z.out.diffx)
#no impact of either treatment

#Trust question in inescapability condition
z.out.trust <- zelig(trust.gov ~ info + threat + bjw + republican, model="ls", data=data)
summary(z.out.trust)
z.out.trustx <- zelig(trust.gov ~ info*threat + bjw + republican, model="ls", data=data)
summary(z.out.trustx)
#no impact of either treatment


